#!/usr/bin/perl -w
use strict;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use uniquify;
use DBI;
use DBD::Pg;
use Data::Dump qw(dump);
use Cwd;
use die_error;
use resolve_gene_to_single_transcript;

#
#	print all discarded ensp if DEBUG == 1
#
use constant DEBUG => 0;






(my $my_dir = $FindBin::Bin) =~ s/(.*)\/$/$1/;



print <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    resolve_ensp_orthologs

	Description:

888888888888888888888888888888888888888888888888888888888888888888888

HEADLINE

#-----------------------------------------------------------------------------------------
#	Get DB parameters
#
use db_parameters;
use ortho_species;
#my @db_param = get_db_parameters();
#
#	Get DB parameters
#-----------------------------------------------------------------------------------------
#
#	initial orthologue transcript predictions (in relative path from here)
#
my $species12 =$ortho_name1."_".$ortho_name2;
my $species1_regex = $ortho_prot1;
my $dir_pipeline_temp_orthologs			= "$dir_pipeline_ds_temp/$species12/orthologs";
my $dir_pipeline_errors					= "$dir_pipeline_ds_temp/$species12/errors";
my $summary_file_name					= "$dir_orthologs_output/resolve_orthologs/resolve.ortholog.counts.summary";
my $sans_loci_error_file_name			= "$dir_pipeline_errors/resolve_orthologs.sequences_sans_loci";
my $ambiguous_file_name					= "$dir_orthologs_output/resolve_orthologs/clades/orthologs_resolve.ambiguous.clades";
my $unambiguous_file_name				= "$dir_orthologs_output/resolve_orthologs/clades/orthologs_resolve.unambiguous.clades";
my $singlechoice_file_name				= "$dir_orthologs_output/resolve_orthologs/clades/orthologs_resolve.single_choice.clades";
my $no_ds_to_ortholog_error_file_name	= "$dir_pipeline_errors/orthology.ortholog_clades.no_ds_to_any_ortholog.errors";
my $ds_pairwise_file_name				= "$dir_orthologs_output/ds_res.orthologs.both_filtered";
###DEBUGGG
#my $ortholog_clades_file_name ="$FindBin::Bin/test.results.initial";
my $ortholog_clades_file_name ="$dir_orthologs_output/interpret_tree/transcript.ortholog.clades";


system ("mkdir -p $dir_pipeline_temp_orthologs");
system ("mkdir -p $dir_pipeline_errors");
system ("mkdir -p $dir_orthologs_output/resolve_orthologs");
system ("mkdir -p $dir_orthologs_output/resolve_orthologs/clades");



#
# connect to genome database
#
my $dbh = connect_to_panda();
$dbh->{RaiseError} = 0;
use get_psql_time;
my $start_time = get_psql_time($dbh);


#_________________________________________________________________________________________
#
#	Get ensp to ensg mapping
#
#_________________________________________________________________________________________
# FUNC_1
sub get_ensp_to_ensg($\%)
{
	my ($dbh, $map_ensp_to_ensg) = @_;
	#ignore previously merged genes
	my $sql_cmd =<<"PL/SQLCMD";
    SELECT
        distinct ON (prot_id)
        gene_id, prot_id
    FROM
        ens_id
    WHERE
        species = '$ortho_name1' OR species = '$ortho_name2'
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd)		or die DBI::errstr;
	print_error("No peptide to gene lookup data retrieved from Panda") unless @$data;
	for (@$data)
	{
		$map_ensp_to_ensg->{$_->[1]} = $_->[0];
	}
}








#_________________________________________________________________________________________
#
#	init_paralog_result_struct
#
#		clade various fields to 0. This is a constructor
#
#_________________________________________________________________________________________
# FUNC_3
sub init_paralog_result_struct(\%$)
{
	my ($paralogs, $name) = @_;

	$paralogs->{NAME}				= $name;

	$paralogs->{COUNT_ORIG_CLADE}	=
	$paralogs->{COUNT_SETS}			= 0;
}






#_________________________________________________________________________________________
#
#	Read in ortholog clades
#
#
#_________________________________________________________________________________________
# FUNC_4
sub read_in_ortholog_clades($\@)
{
	my ($file, $all_clades) = @_;
	open_or_die (*FILE, $file);
	for (<FILE>)
	{
		chomp;
		#ignore blank or header lines
		next if (!length($_) or /^>/);

		my ($dS,  @ensp) = split /\t/;

		#save orthologue clade
		push(@$all_clades, [@ensp]);

		if (DEBUG)
		{
			print "clade_id = ", ( @$all_clades - 1), "\t", join (", ",  sort @ensp), "\n";
		}
	}
}



#_________________________________________________________________________________________
#
#	split_transcripts_by_species
#
#_________________________________________________________________________________________
# FUNC_5
sub split_transcripts_by_species(\@)
{
	my ($all_clades) = @_;
	for my $ensp(@$all_clades)
	{
		#
		# divide ensp per species
		#
		my (@ensp1,  @ensp2);
		for (@$ensp)
		{
			if (/$species1_regex/)
			{
				push(@ensp1,  $_);
			}
			else
			{
				push(@ensp2,  $_);
			}
		}
		$ensp = [[@ensp1], [@ensp2]];
	}
}








#_________________________________________________________________________________________
#
#	record_unclustered_clades
#
#_________________________________________________________________________________________
# FUNC_6
sub record_unclustered_clades(\@\%)
{
	my ($all_clades, $results) = @_;
	for my $clade ( @$all_clades)
	{
		my ($ensp1, $ensp2) = @$clade;
		# save transcripts
		@{$results->{0}{ORIG_TRANSCRIPTS}}{@$ensp1} = ();
		@{$results->{1}{ORIG_TRANSCRIPTS}}{@$ensp2} = ();
	}
	$results->{COUNT_ORIG_CLADE}	=
	$results->{COUNT_SETS}			= scalar @$all_clades;
}









#_________________________________________________________________________________________
#
#	remove_overlapping_transcripts
#
#		removing all overlapping transcripts
#
#_________________________________________________________________________________________
# FUNC_7.21
my %ensp_sans_loci;
sub remove_overlapping_transcripts($\%\%\%\%\%)
{
	my ($ensp, 
		$ensg_to_ensp, 
		$ensp_to_ensg,
		$ensp_to_loci, 
		$discarded_ids,
		$unresolved_ids) = @_;

	# compare this transcript with all others from the same gene
	my $ensg = $ensp_to_ensg->{$ensp};
	my @ids = keys %{$ensg_to_ensp->{$ensg}};
	for my $id(@ids)
	{
		# avoid self-self comparisons
		next if $id eq $ensp;

		# delete if overlap
		# FUNC_7.211
		if (overlap($id, $ensp, %$ensp_to_loci, %ensp_sans_loci))
		{
			++$discarded_ids->{$id};
			delete $unresolved_ids->{$id};
			delete $ensg_to_ensp->{$ensg}{$id};
		}
	}
	delete $unresolved_ids->{$ensp};
}

#_________________________________________________________________________________________
#
#	keep_non_overlapping_transcripts_with_lowest_ds
#
#		1. start from lowest dS pairs
#		2. delete all overlapping for each identifier of pair (discarded_ids)
#		3. use next lowest dS
#
#		Overlapping transcripts
#			1. Removed from $ensg_to_ensp
#			2. Added to $discarded_ids
#			3. removed from $unresolved_idsXX
#		Kept transcripts
#			1. Removed from $unresolved_idsXX
#			2. Added to $resolved_idsXX
#
#_________________________________________________________________________________________
# FUNC_7.2
use constant ID1 => 1;
use constant ID2 => 2;
sub keep_non_overlapping_transcripts_with_lowest_ds(\@\%\%\%\%\%\%\%\%\%)
{
	my (
			$ds,				# dS pairwise relationships [dS, id1, id2]
								# id1 and id2 are either	species 1 / species 2 or 
								#							resolved  / unresolved
			$discarded_ids1,	# overlapping transcripts with greater dS
			$discarded_ids2,	# overlapping transcripts with greater dS
			$ensg_to_ensp, 		# map ensg to array of ensp
			$ensp_to_ensg,		# map ensp back to ensg
			$ensp_to_loci,		# map ensp to its loci
			$unresolved_ids1,	# list of unresolved ids
			$unresolved_ids2,
			$resolved_ids1,		# list of resolved ids
			$resolved_ids2,

			) = @_;

	my $cnt_resolved = 0;
	for my $i (0.. (@$ds - 1))
	{
		my ($ensp1, $ensp2) = ($ds->[$i][ID1], $ds->[$i][ID2]);

		# ignore if the relationship is with a discarded transcript
		next if exists $discarded_ids1->{$ensp1} ||
				exists $discarded_ids2->{$ensp2};
	
		# find and discard all overlapping for transcripts 1
		if (exists $unresolved_ids1->{$ensp1})
		{
			# FUNC_7.21
			remove_overlapping_transcripts($ensp1, %$ensg_to_ensp, %$ensp_to_ensg,
										   %$ensp_to_loci, 
										   %$discarded_ids1, 
										   %$unresolved_ids1);
			++$resolved_ids1->{$ensp1};
			++$cnt_resolved;
		}
	
		# find and discard all overlapping for transcripts 2
		if (exists $unresolved_ids2->{$ensp2})
		{
			# FUNC_5.21
			remove_overlapping_transcripts($ensp2, %$ensg_to_ensp, %$ensp_to_ensg,
										   %$ensp_to_loci, 
										   %$discarded_ids2, 
										   %$unresolved_ids2);
			++$resolved_ids2->{$ensp2};
			++$cnt_resolved;
		}
	}
	return 	$cnt_resolved;
}




#_________________________________________________________________________________________
#
#	resolve_overlapping_transcripts_within_all_clades
#
#		each orthologue clade must not contain overlapping transcripts from multiple
#			genes
#		resolve by taking transcript with lowest ortholog dS, removing all overlapping
#			and recurse
#
#_________________________________________________________________________________________

#	save 
#			# of starting transcripts
#			# of transcripts with orthology relationships
#			# of transcripts without orthology relationships
#			# resolved/unresolved of each
#
#	DISCARDED_OVERLAPPING
#	DISCARDED_ORPHANS
#	VIA_PARALOGS
#
# FUNC_7
sub resolve_overlapping_transcripts_within_all_clades(\@\%\%\%\%)
{
	my ($clades, 
		$ensp_to_ensg, 
		$ensp_to_loci,
		$seq_pairs_dS,
		$ortholog_results) = @_;

	for my $clade(@$clades)
	{
		#
		# ensp per species
		#
		my ($all_ensp1,  $all_ensp2) = @$clade;

		#
		#   allow ensg to ensp lookup
		#
		my (%ensg_to_ensp);
		for my $ensp(@$all_ensp1, @$all_ensp2)
		{
			die $ensp unless exists $ensp_to_ensg->{$ensp};
			
			my $ensg = $ensp_to_ensg->{$ensp};
			++$ensg_to_ensp{$ensg}{$ensp};
		}

		# overlapping identifiers with larger dS
		my @discarded_ids;
		# identifiers to keep
		my @resolved_ids;
		# identifiers to be resolved
		my @unresolved_ids;
		++$unresolved_ids[0]{$_} for (@$all_ensp1);
		++$unresolved_ids[1]{$_} for (@$all_ensp2);

		#
		#	Deal with orthologous relationships
		#		
		#		1.	keep identifiers starting from lowest dS
		#		2.	remove overlapping transcripts
		{
			# FUNC 7.1
			# get dS relationship
			my @ds = get_sorted_pairswise_ds(@$all_ensp1, @$all_ensp2, %$seq_pairs_dS);
			# FUNC 7.2
			keep_non_overlapping_transcripts_with_lowest_ds(@ds, 
															%{$discarded_ids[0]},
															%{$discarded_ids[1]},
															%ensg_to_ensp,
															%$ensp_to_ensg,
															%$ensp_to_loci,
															%{$unresolved_ids[0]},
															%{$unresolved_ids[1]},
															%{$resolved_ids[0]},
															%{$resolved_ids[1]});

			#
			#	save transcripts with orthologous relationships
			#
			$clade = [[keys %{$resolved_ids[0]}], [keys %{$resolved_ids[1]}]];
		}					
		#
		#	Deal with transcripts without orthologous relationships
		#		
		#		1.	keep identifiers starting from lowest dS
		#		2.	remove overlapping transcripts
		for my $species(0,1)
		{
			while (1)
			{
				#
				# get list of dS paralogous relationships between resolved and unresolved sorted by dS
				#
				my @all_resolved_ids = keys %{$resolved_ids[$species]};
				# start resolved_ids from scratch so can save paralogs
				$resolved_ids[$species] = {};
				my @all_unresolved_ids = keys %{$unresolved_ids[$species]};
				my @ds = get_sorted_pairswise_ds(	@all_resolved_ids, 
													@all_unresolved_ids, 
													%$seq_pairs_dS);


				my $cnt_resolved = 
					keep_non_overlapping_transcripts_with_lowest_ds(@ds, 
																	%{$discarded_ids[$species]},
																	%{$discarded_ids[$species]},
																	%ensg_to_ensp,
																	%$ensp_to_ensg,
																	%$ensp_to_loci,
																	%{$unresolved_ids[$species]},
																	%{$unresolved_ids[$species]},
																	%{$resolved_ids[$species]},
																	%{$resolved_ids[$species]});
				# break if nothing more to resolve (leaving orphans) 
				# otherwise recurse
				last unless $cnt_resolved;
				last unless scalar %{$unresolved_ids[$species]};
			}

			#
			#	save all transcripts with paralogous relationships for this species
			#
			for (keys %{$resolved_ids[$species]})
			{
				push(@{$clade->[$species]}, $_);
			}
		}



		for my $species(0,1)
		{
			# save discarded identifiers representing overlapping transcripts
			for (keys %{$discarded_ids[$species]})
			{
				++$ortholog_results->{$species}{DISCARDED_OVERLAPPING}{$_};
			}

			# save discarded identifiers representing orphaned transcripts
			for (keys %{$unresolved_ids[$species]})
			{
				++$ortholog_results->{$species}{DISCARDED_ORPHANS}{$_};
			}

			# save discarded identifiers representing orphaned transcripts
			for (keys %{$resolved_ids[$species]})
			{
				++$ortholog_results->{$species}{VIA_PARALOGS}{$_};
			}
		}


		
	}

}





#_________________________________________________________________________________________


#	calculate_mean_ortholog_dS_for_all_clades

#		for an ortholog clades (of ensp), get the average dS for orthology comparisons


#_________________________________________________________________________________________
# FUNC_8.1
sub calculate_mean_ortholog_dS(\@\@\%)
{
	my ($species1_ensp, $species2_ensp,  $seq_pairs_dS) = @_;
	# all vs all dS
	my $cnt_dS = 0;
	my $total_dS = 0;
	for my $ensp1(@$species1_ensp)
	{
		for my $ensp2(@$species2_ensp)
		{
			if (exists $seq_pairs_dS->{$ensp1} &&
				exists $seq_pairs_dS->{$ensp1}{$ensp2})
			{
				$total_dS += $seq_pairs_dS->{$ensp1}{$ensp2};
				++$cnt_dS;
			}
		}
	}
	# fail: no dS
	if ($cnt_dS == 0)
	{
		return 10;
	}
	return $total_dS / $cnt_dS;
}
# FUNC_8
sub calculate_mean_ortholog_dS_for_all_clades(\@\%\@)
{
	my ($all_clades, $seq_pairs_dS, $mean_ortholog_dS) = @_;
	for my $ensp(@$all_clades)
	{
		push(@$mean_ortholog_dS, calculate_mean_ortholog_dS(@{$ensp->[0]},
															@{$ensp->[1]},
															%$seq_pairs_dS));  
	}
}



















#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	assign_ensp_to_clade_with_lowest_ds
#		also calculate mean ortholog ds and splits transcripts by species
#
#
# FUNC_9.1
sub delete_item_from_array(\@$)
{
	my ($arr, $item) = @_;
	my %items;
	@items{@$arr}= ();
	delete $items{$item};
	@$arr = keys %items;
	
}
# FUNC_9
sub assign_ensp_to_clade_with_lowest_ds(\@\@)
{
	my ($mean_ortholog_ds, $all_clades) = @_;
	
	my %ensp_to_clade_id;
	my %ensp_to_dS;

	#
	#	For each transcript in each clade, 
	#		test if already in another clade?
	#			Yes: Does the transcript have a lower dS there?
	#				Yes: Discard the transcript in this clade
	#			No: Discard the transcript in the other clade
	#			
	#
	for my $clade_id(0.. (@$all_clades - 1))
	{
		my $clade_ds = $mean_ortholog_ds->[$clade_id];

		my $clade = $all_clades->[$clade_id];
		for my $species_id(0,1)
		{
			my $ensp_in_clade = $clade->[$species_id];

			#only keep ensp with lowest dS
			my @ensp_to_keep_in_this_clade;

			for my $ensp(@$ensp_in_clade)
			{
				if (exists $ensp_to_dS{$ensp})
				{
					# keep the saved ensp in the other clade if it has a lower dS
					next if ($clade_ds > $ensp_to_dS{$ensp});

					# otherwise delete the saved ensp and keep this
					my $other_clade_id = $ensp_to_clade_id{$ensp};
					# FUNC_9.2
					delete_item_from_array(@{$all_clades->[$other_clade_id][$species_id]}, $ensp);
				}

				# save this ensp for this clade
				$ensp_to_dS{$ensp} = $clade_ds;
				$ensp_to_clade_id{$ensp} = $clade_id;
				push(@ensp_to_keep_in_this_clade, $ensp);
			}
			@{$clade->[$species_id]} = @ensp_to_keep_in_this_clade;
		}
	}
}





#_________________________________________________________________________________________
#
#	cluster_ortholog_clades(\@$\%) 
#
#		cluster ortholog_clades which share common ensg
#
#_________________________________________________________________________________________
# FUNC_10
sub cluster_ortholog_clades(\@\@$\%)
{
	my ($clusters, $unlinked_clades, $size, $ensg_to_clade_id) = @_;
	my @clusters_of_clades;
	my %clade_id_to_cluster_id;

	# start with each orthoclade in its own cluster
	for (0..($size - 1))
	{
		push(@clusters_of_clades, [$_]);
		$clade_id_to_cluster_id{$_} = $_;
	}
	

	# for each ensg
	for my $ensg(keys %$ensg_to_clade_id)
	{
		# get all the orthoclades it is a member of and link them
		my $clade_ids = $ensg_to_clade_id->{$ensg};

		# this gene only found in one orthoclade
		next if @$clade_ids == 1;


		# the cluster of the first orthoclade is where everything will end up
		my $cluster_id0 = $clade_id_to_cluster_id{$clade_ids->[0]};
		for my $clade_id(@$clade_ids)
		{
			my $orig_cluster_id = $clade_id_to_cluster_id{$clade_id};
			next if ($orig_cluster_id == $cluster_id0);

			#
			# merge contents of $orig_cluster_id with $cluster_id0
			#

			# 1. old clade_ids should point to cluster_id0
			for my $orig_clade_id (@{$clusters_of_clades[$orig_cluster_id]})
			{
				$clade_id_to_cluster_id{$orig_clade_id} = $cluster_id0;
			}

			# 2. move contents of $orig_cluster_id to $cluster_id0;
			push(@{$clusters_of_clades[$cluster_id0]}, @{$clusters_of_clades[$orig_cluster_id]});
			@{$clusters_of_clades[$orig_cluster_id]} = ();
		}
	}

	# only put non empty clusters into results
	for (@clusters_of_clades)
	{
		if (@$_ > 1)
		{
			push(@$clusters, $_) 
		}
		elsif (@$_ == 1)
		{
			push(@$unlinked_clades, @$_) 
		}
	}
}







#_________________________________________________________________________________________


#	ortholog_clades_in_cluster_are_identical

#		for a cluster of clades, check if all ortholog clades are identical


#_________________________________________________________________________________________
# FUNC_11
sub ortholog_clades_in_cluster_are_identical(\@\@)
{
	my ($cluster_of_clade_ids, $ortholog_clades) = @_;
	return 1 if (@$cluster_of_clade_ids == 1);
	#
	# get string representation of ortholog clades in cluster
	# 
	my @ortholog_clade_str;
	for my $clade_id(@$cluster_of_clade_ids)
	{
		push(@ortholog_clade_str, join("??",	@{$ortholog_clades->[$clade_id][0]},
												@{$ortholog_clades->[$clade_id][1]}));
	}

	# compare first ortholog clade with the rest
	for my $i (1 .. $#ortholog_clade_str)
	{
		return 0 if $ortholog_clade_str[0] ne $ortholog_clade_str[$i];
	}
	return 1;
}








#_________________________________________________________________________________________
#
#	calculate_lowest_ortholog_dS_per_ensp
#
#		For each ensp, calculate the lowest dS to any orthologue
#
#_________________________________________________________________________________________
# FUNC_12.111
sub get_ortholog_ds_from_closest_paralog(\@\%\%)
{
	my ($all_ensp, $seq_pairs_dS, $ortholog_ds) = @_;
	# for those without, use closest paralog
	TRANSCRIPT:
	for my $ensp1(@$all_ensp)
	{
		next if (exists $ortholog_ds->{$ensp1});

		# get all paralogs for this transcript
		my @ds;
		for my $ensp2(@$all_ensp)
		{
			next if $ensp1 eq $ensp2;
			next unless exists $seq_pairs_dS->{$ensp1}{$ensp2};
			my $ds = $seq_pairs_dS->{$ensp1}{$ensp2};
			push(@ds,  [$ensp2, $ds]);
		}

		# sort by ds
		@ds = sort {$a->[1] <=> $b->[1]} @ds;

		# go through paralogs and use their orthologous distance if any
		for my $ensp_ds(@ds)
		{
			my ($ensp, $ds) = @$ensp_ds;
			if (exists $ortholog_ds->{$ensp})
			{
				$ortholog_ds->{$ensp1} = $ortholog_ds->{$ensp};
				next TRANSCRIPT; 
			}
		}
	}
}
# FUNC_12.11
sub calculate_lowest_ortholog_dS_per_ensp(\@\@\%\%)
{
	my ($species1_ensp, $species2_ensp, $seq_pairs_dS, $ortholog_ds) = @_;

	# get lowest ortholog ds
	for my $ensp1(@$species1_ensp)
	{
		for my $ensp2(@$species2_ensp)
		{
			if (exists $seq_pairs_dS->{$ensp1}{$ensp2})
			{
				my $ds = $seq_pairs_dS->{$ensp1}{$ensp2};
				if (!exists $ortholog_ds->{$ensp1})
				{
					$ortholog_ds->{$ensp1} = $ds;
				}
				else
				{
					$ortholog_ds->{$ensp1} = $ds if ($ds < $ortholog_ds->{$ensp1});
				}
				if (!exists $ortholog_ds->{$ensp2})
				{
					$ortholog_ds->{$ensp2} = $ds;
				}
				else
				{
					$ortholog_ds->{$ensp2} = $ds if ($ds < $ortholog_ds->{$ensp2});
				}
			}
		}
	}
	# if does not have lowest ortholog ds
	#		use one from a paralog
	# FUNC_12.111
	get_ortholog_ds_from_closest_paralog(@$species1_ensp, %$seq_pairs_dS, %$ortholog_ds);
	get_ortholog_ds_from_closest_paralog(@$species2_ensp, %$seq_pairs_dS, %$ortholog_ds);
}


#_________________________________________________________________________________________
#
#	print_and_record_clade
#
#_________________________________________________________________________________________
# FUNC_12.1
my $ortholog_id = 0;
sub validate_and_print_clade(\@*\%)
{
	my ($transcripts,
		$FH_RESULTS,
		$seq_pairs_ds) = @_;



	#
	#	print out results
	#
	# 	FUNC_12.11
	my %lowest_ortholog_ds;
	calculate_lowest_ortholog_dS_per_ensp(@{$transcripts->[0]}, @{$transcripts->[1]}, %$seq_pairs_ds, 
										%lowest_ortholog_ds);
	my @validated_transcripts;
	for my $species(0,1)
	{
		for my $transcript (@{$transcripts->[$species]})
		{
			next unless exists $lowest_ortholog_ds{$transcript};
			print $FH_RESULTS join ("\t",	$ortholog_id, 
											$lowest_ortholog_ds{$transcript},
											$transcript), "\n";
			push(@{$validated_transcripts[$species]}, $transcript);
		}
	}
	++$ortholog_id;

	$transcripts = [@validated_transcripts];
}


sub store_clade_transcripts(\@\%$)
{
	my ($transcripts, 
		$ortholog_results,
		$category) = @_;

	#
	#	count genes and paralogs
	#
	for my $species(0,1)
	{
		# save transcripts
		@{$ortholog_results->{$species}{$category}}{@{$transcripts->[$species]}} = ();
	}
}



#_________________________________________________________________________________________
#
#	resolve_singly_spliced_transcripts
#
#		Add annotation for merged genes in panda
#
#_________________________________________________________________________________________
# FUNC_12
sub resolve_singly_spliced_transcripts(\@\@\@\%\%\%)
{
	my ($unlinked_ortho_clades, 
		$mean_ortholog_ds,
		$all_clades,
		$map_ensp_to_ensg,
		$seq_pairs_dS,
		$ortholog_results) = @_;

	open_or_die (*FH_RESOLVED, ">$singlechoice_file_name");


	for my $clade_id (@$unlinked_ortho_clades)
	{
		#
		#	error if no orthologous links
		#
		if ($mean_ortholog_ds->[$clade_id] == 10)
		{
			# should not happen because filtered already
			# Bug
			print STDERR "Bug. No orthologue dS between any members!\t".
							"clade_id = $clade_id\tMembers of clade=>\t",
							join (", ",	@{$all_clades->[$clade_id][0]}, 
										@{$all_clades->[$clade_id][1]}), "\n";
			die if ($mean_ortholog_ds->[$clade_id] == 10);
		}

		validate_and_print_clade(	@{$all_clades->[$clade_id]},
									*FH_RESOLVED,
									%$seq_pairs_dS);
		store_clade_transcripts(	@{$all_clades->[$clade_id]},
									%$ortholog_results,
									'TRANSCRIPTS');
	}
	$ortholog_results->{0}{'ORIG_TRANSCRIPTS'} = $ortholog_results->{0}{'TRANSCRIPTS'};
	$ortholog_results->{1}{'ORIG_TRANSCRIPTS'} = $ortholog_results->{1}{'TRANSCRIPTS'};
	
}








#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	count_clades_in_sets
#
#
# FUNC_13
sub count_clades_in_sets(\%\@)
{
	my ($ortholog_results, $clade_sets) = @_;
	$ortholog_results->{COUNT_ORIG_CLADE}	= 0;
	for my $clade_set(@$clade_sets)
	{
		$ortholog_results->{COUNT_ORIG_CLADE} += @$clade_set;
	}
	$ortholog_results->{COUNT_SETS} = scalar @$clade_sets;
}







#_________________________________________________________________________________________


#	ortholog_clade_with_lowest_ds

#		for a cluster of clades, get the ortholog clade with the lowest dS


#_________________________________________________________________________________________
# FUNC_14.2
sub ortholog_clade_with_lowest_ds(\@\@)
{
	my ($cluster_of_clade_ids, $all_clades_mean_ortholog_ds) = @_;

	my $lowest_dS = 9999.0;
	my $best_clade_id;
	for my $clade_id(@$cluster_of_clade_ids)
	{
		die join ("\n", @$cluster_of_clade_ids) if $clade_id eq '';
		my $mean_ortholog_ds = $all_clades_mean_ortholog_ds->[$clade_id];
		
		if ($mean_ortholog_ds < $lowest_dS)
		{
			$lowest_dS = $mean_ortholog_ds;
			$best_clade_id = $clade_id;
		}
	}
	return $best_clade_id;
}







#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	record_paralog_count_per_clade_set
#
#
# FUNC_14.1
sub record_paralog_count_per_clade_set($$$$)
{
	my ($clade_set,  
		$all_clades, 
		$map_ensp_to_ensg,
		$paralogs) = @_;


	for my $species(0,1)
	{
		# save all transcripts
		for my $clade_id(@$clade_set)
		{
			my $ensp = $all_clades->[$clade_id][$species];
			@{$paralogs->{$species}{ORIG_TRANSCRIPTS}}{@$ensp} = ();
		}
	}
}


#_________________________________________________________________________________________
#
#	resolve_unambiguous_clades
#
#		resolve transcripts clades which all have corresponding genes and vice versa
#
#_________________________________________________________________________________________
# FUNC_14
sub resolve_unambiguous_clades($\@\@\%*\%\%)
{
	my ($clade_ids,  
		$mean_ortholog_ds, 
		$all_clades, 
		$map_ensp_to_ensg,
		$FH_RESOLVED,
		$seq_pairs_dS,
		$ortholog_results,
		) = @_;


	# FUNC_14.1
	record_paralog_count_per_clade_set($clade_ids, $all_clades, $map_ensp_to_ensg,
										$ortholog_results);
	
	# FUNC_14.2
	# will blow up with undef if this clade_id has no mean dS
	my $clade_id = ortholog_clade_with_lowest_ds(@$clade_ids, @$mean_ortholog_ds);

	#
	#	error if no orthologous links
	#
	die if ($mean_ortholog_ds->[$clade_id] == 10);
	if ($mean_ortholog_ds->[$clade_id] == 10)
	{
		# should not happen because filtered already
		# Bug
		die "Error:\n\tBug. No orthologue dS between any members!\t".
						"clade_id = $clade_id\tMembers of clade=>\t",
						join (", ",	@{$all_clades->[$clade_id][0]}, 
									@{$all_clades->[$clade_id][1]}), "\n";
	}

	validate_and_print_clade(	@{$all_clades->[$clade_id]},
								$FH_RESOLVED,
								%$seq_pairs_dS);
	store_clade_transcripts(	@{$all_clades->[$clade_id]},
								%$ortholog_results,
								'TRANSCRIPTS');

	if (DEBUG)
	{
		my $best_clade_id = $clade_id;
		for my $clade_id (@$clade_ids)
		{
			next if $clade_id == $best_clade_id;
			print "Unambiguous clade ($clade_id) discarded (in favour of $best_clade_id)\t:", join (", ",
												   @{$all_clades->[$clade_id][0]}, 
                                                   @{$all_clades->[$clade_id][1]}), 
				  "\n";
		}
	}
}






#_________________________________________________________________________________________
#
#	resolve_ambiguous_clades
#
#
#		Break apart (Resolve) a clade of orthologue clades which have overlapping and
#			different number of transcripts in different clades. 
#		There are probably merged genes as well
#
#		For ortholog clades which have been linked together because they share ensg
#			1) keep the clade which has the lowest ortholog dS
#			2) Remove all other clades which
#						a) have shared ensp
#						b) have any ensp with the same gene which overlap any from the
#								lowest dS clade
#			3) Repeat
#
#_________________________________________________________________________________________
# FUNC_15
sub resolve_ambiguous_clades($\@\@*\%\%\%\%)
{
	my ($clades,  
		$mean_ortholog_ds, 
		$all_clades, 
		$FH_RESOLVED,
		$map_ensp_to_ensg,
		$map_ensp_to_loci,
		$seq_pairs_dS,
		$ortholog_results) = @_;

	record_paralog_count_per_clade_set($clades, $all_clades, $map_ensp_to_ensg,
										$ortholog_results);

	#
	# 	Sort ensp by ensg 
	#		So that alternative splicing / presence in multiple clades can
	#			be easily detected. 
	#		Use hash rather than array for easier deletes
	#
	my @clade_genes_to_transcripts;
	my @mean_ortholog_ds;
	my $cnt_ortholog_clades = 0;

	my $logging = 0;
	
	for my $clade_id(@$clades)
	{
		###DEBUGGG
		#$logging = 1 if ($clade_id == 1242);

		# save starting dS
		push(@mean_ortholog_ds, $mean_ortholog_ds->[$clade_id]);

		for my $species(0,1)
		{
			for my $ensp(@{$all_clades->[$clade_id][$species]})
			{
				my $ensg = $map_ensp_to_ensg->{$ensp};
				$clade_genes_to_transcripts[$cnt_ortholog_clades][$species]{$ensg}{$ensp}++;
			}
		}
		++$cnt_ortholog_clades;
	}
	print dump @clade_genes_to_transcripts if $logging;
	

	# start off with all orthologue clades
	my %unresolved_clades;
	++$unresolved_clades{$_} for (0..($cnt_ortholog_clades - 1));
	

	# new clusters
	my @resolved_clade_data;
	my @resolved_clades;

	my $iteration = 0;
	
	# 
	# Resolving
	#
	while (1)
	{
		# stop if all resolved
		last unless keys %unresolved_clades;

		###LOGGING
		++$iteration;
		print "\nITERATION $iteration\n" if $logging;

		#
		#	Get unresolved ortholog clade with lowest dS
		#
		my @unresolved_clades = keys %unresolved_clades;
		my $best_clade_id = ortholog_clade_with_lowest_ds(@unresolved_clades, 
														   @mean_ortholog_ds);
		my $best_ortholog_clade = $clade_genes_to_transcripts[$best_clade_id];

		###LOGGING
		print "Best ortholog clade = \n", dump $best_ortholog_clade if $logging;
		
		#
		#	This clade now 'resolved'
		#
		push(@resolved_clades, $best_ortholog_clade);
		delete($unresolved_clades{$best_clade_id});

		###LOGGING
		print "Resolved clades = \n", dump @resolved_clades if $logging;
		print "Unresolved clades = \n", dump %unresolved_clades if $logging;

		
		#  and save mean ortholog dS and ensp
		my @ensp;
		@{$ensp[0]} = ();
		@{$ensp[1]} = ();
		for my $species(0, 1)
		{
			for my $ensg (keys %{$best_ortholog_clade->[$species]})
			{
				push(@{$ensp[$species]}, keys %{$best_ortholog_clade->[$species]{$ensg}});
			}
		}
		# both species must have transcripts!
		die (join (",", @$clades)) unless scalar @{$ensp[0]} && scalar @{$ensp[1]};
		#
		#	error if no orthologous links
		#
		die "$best_clade_id:", (join (",", @$clades, @{$ensp[0]}, @{$ensp[1]}, @mean_ortholog_ds)) 
									if $mean_ortholog_ds[$best_clade_id] >= 5;
		push(@resolved_clade_data, [@ensp]);

		###LOGGING
		print "Resolved clade data = \n", dump @resolved_clade_data if $logging;


		# for each of the remaining ortholog clades
		# remove all overlapping transcripts
		# remove if ds > 5
		# remove if not both species populated

		@unresolved_clades = keys %unresolved_clades;
		SET_LOOP:
		for my $clade_id(@unresolved_clades)
		{	
			my $clade_modified = 0;

			# for each gene
			for my $species(0, 1)
			{
				my @shared_ensp;
				my @overlapping_ensp;
				
				# an entire clade of genes for this species
				my $ortholog_clade = $clade_genes_to_transcripts[$clade_id][$species];

				for my $ensg (keys %{$best_ortholog_clade->[$species]})
				{
					# Is this gene cross linking this clade?
					#	i.e. Does this gene exists in this clade as well?
					next unless exists $ortholog_clade->{$ensg};

					# get peptides for this gene in this clade 
					my $ortholog_clade_ensp = $ortholog_clade->{$ensg};
	
					
					# Remove all transcripts which are shared with or overlap those
					#	from the "best" clade
					# 	The old behaviour was to delete the entire clade
					my @ensp_resolved	= keys %{$best_ortholog_clade->[$species]{$ensg}};
					my @ensp			= keys %$ortholog_clade_ensp;
					for my $ensp_resolved (@ensp_resolved)
					{
						for my $ensp(@ensp)
						{
							if ($ensp eq $ensp_resolved)
							{
								###LOGGING
								print "Matching transcripts with best clade $ensp in clade=",
										dump $ortholog_clade, "\n" if $logging;

								push(@shared_ensp, $ensp);
							}
							elsif (overlap($ensp, $ensp_resolved, %$map_ensp_to_loci, %ensp_sans_loci))
							{
								###LOGGING
								print "Overlapping transcript [$ensp] with best clade [$ensp_resolved] in clade=",
										dump $ortholog_clade, "\n" if $logging;

								push(@overlapping_ensp, $ensp);
							}
							else
							{
								next;
							}

							###LOGGING
							print "Delete $ensp \n" if $logging;

									# if shared or overlap: delete
							delete $ortholog_clade_ensp->{$ensp};
							++$clade_modified;

							# old behaviour:delete entire clade if overlap
							#   delete($unresolved_clades{$clade_id});
							#   next SET_LOOP;
						}
	
					}

					if (DEBUG || $logging)
					{
						if (@overlapping_ensp)
						{
							print "Ambiguous clade id = $clade_id: overlapping transcripts vs ".
									"best id = $best_clade_id for species $species\t",
										join (", ", @overlapping_ensp), "\n";
						}
						if (@shared_ensp)
						{
							print "Ambiguous clade id = $clade_id: shared transcripts vs ".
									"best id = $best_clade_id for species $species\t",
										join (", ", @shared_ensp), "\n";
						}
					}
	
					# if all transcript have been removed from this clade, remove gene as well
					if (!scalar keys %$ortholog_clade_ensp)
					{
						###LOGGING
						print "No more transcripts from $ensg remove Gene\n" if $logging;

						delete $ortholog_clade->{$ensg};
						#delete $ortholog_clade_ensp;
					}
				}
			}

			#
			#	delete clade if does not have genes/transcripts from both species
			#
			my $ortholog_clade = $clade_genes_to_transcripts[$clade_id];
			my $ortholog_clade1 = $ortholog_clade->[0];
			my $ortholog_clade2 = $ortholog_clade->[1];
			if (keys %$ortholog_clade1 == 0 ||
				keys %$ortholog_clade2 == 0)
			{
				if (DEBUG || $logging )
				{
					print "Ambiguous clade discard clade $clade_id with only ".
							"single species\t", 
					join(", ", keys %$ortholog_clade1, keys %$ortholog_clade2), "\n";
				}

				delete($unresolved_clades{$clade_id});
				next;
			}

			#
			#	recalculate mean dS if clade modified
			#
			if ($clade_modified)
			{
				# obtain all ensp in the clade
				my @ensp;
				for my $species(0, 1)
				{
					for my $ensg (keys %{$ortholog_clade->[$species]})
					{
						push(@{$ensp[$species]},  keys %{$ortholog_clade->[$species]{$ensg}});
					}
				}

				$mean_ortholog_ds[$clade_id] = 
					calculate_mean_ortholog_dS(@{$ensp[0]}, @{$ensp[1]}, %$seq_pairs_dS);
				###LOGGING
				print "mean_ortholog_ds[$clade_id]=",$mean_ortholog_ds[$clade_id],
					join (", ",  @{$ensp[0]}, @{$ensp[1]}), "\n" if $logging;
			}

			# only keep clades with ortholog dS <= 5
			if ($mean_ortholog_ds[$clade_id] >= 5)
			{
				if (DEBUG || $logging)
				{
					print "Ambiguous clade discard high dS clade $clade_id\t", 
					join(", ", keys %$ortholog_clade1, keys %$ortholog_clade2), "\n";
				}
				delete($unresolved_clades{$clade_id});
				next;
			}
		}	# for each ortholog clade
		
	}	# recurse until all resolved


	print "Final\n", dump @resolved_clade_data, "\n" if $logging;
	#
	# print out ensp ortholog clades
	#
	for my $i (0 .. $#resolved_clade_data)
	{	
		validate_and_print_clade(	@{$resolved_clade_data[$i]},
									$FH_RESOLVED,
									%$seq_pairs_dS);
		store_clade_transcripts(	@{$resolved_clade_data[$i]},
									%$ortholog_results,
									'TRANSCRIPTS');

		# store additional clades as 'extra' transcripts
		next if $i;
		store_clade_transcripts(	@{$resolved_clade_data[$i]},
									%$ortholog_results,
									'NON_EXTRA_TRANSCRIPTS');
	}

}
















#_________________________________________________________________________________________
#
#	merge_paralog_result_struct
#
#		merge the results from paralogs2 to paralogs1
#
#_________________________________________________________________________________________
# FUNC_16
sub merge_paralog_result_struct(\%\%)
{
	my ($paralogs1, $paralogs2) = @_;
	# add counts
	for my $category(qw(COUNT_SETS COUNT_ORIG_CLADE) )
	{
		$paralogs1->{$category}	+= $paralogs2->{$category}; 
	}
	for my $species(0,1)
	{
		# add transcripts
		for my $category (qw(TRANSCRIPTS ORIG_TRANSCRIPTS))
		{
			my $paralogs1_sp = \%{$paralogs1->{$species}{$category}};
			my $paralogs2_sp = \%{$paralogs2->{$species}{$category}};
			@{$paralogs1_sp}{keys %$paralogs2_sp} = ();
		}
	}
}




#_________________________________________________________________________________________
#
#	finish_transcript_gene_counts
#
#_________________________________________________________________________________________
# FUNC_17
sub finish_transcript_gene_counts(\%\%)
{
	my ($counts_all_species, $map_ensp_to_ensg) = @_;
	for my $species(0,1)
	{
		my $counts = $counts_all_species->{$species};

		#
		#	get genes from transcripts
		#
		$counts->{GENES}{$map_ensp_to_ensg->{$_}}{$_} = ()
								for (keys %{$counts->{TRANSCRIPTS}});
		$counts->{ORIG_GENES}{$map_ensp_to_ensg->{$_}}{$_} = ()
								for (keys %{$counts->{ORIG_TRANSCRIPTS}});


		#
		#	initialize counts to zero
		#
		$counts->{COUNT_GENES_MULTIPLE_TRANSCRIPTS}	=
		$counts->{COUNT_MULTIPLE_TRANSCRIPTS}		=
		$counts->{COUNT_GENES_EXTRA_TRANSCRIPTS}	= 0;
			

		#
		#	Do counts
		#
		
		# count genes and transcripts
		$counts->{COUNT_GENES}				= scalar keys %{$counts->{GENES}};
		$counts->{COUNT_TRANSCRIPTS}		= scalar keys %{$counts->{TRANSCRIPTS}};
		$counts->{COUNT_ORIG_GENES}			= scalar keys %{$counts->{ORIG_GENES}};
		$counts->{COUNT_ORIG_TRANSCRIPTS}	= scalar keys %{$counts->{ORIG_TRANSCRIPTS}};
	
		# count the number of genes ambiguous would have predicted without demerging
		my %non_extra_genes	= map {$map_ensp_to_ensg->{$_} => undef} 
										keys %{$counts->{NON_EXTRA_TRANSCRIPTS}};
		my $cnt_non_extra_genes = scalar keys %non_extra_genes;

		$counts->{COUNT_EXTRA_TRANSCRIPTS} = $cnt_non_extra_genes ?
							$counts->{COUNT_TRANSCRIPTS} - $cnt_non_extra_genes :
							$counts->{COUNT_TRANSCRIPTS} - $counts->{COUNT_GENES};
							


		# count genes with multiple transcripts
		for my $gene_id(keys %{$counts->{GENES}})
		{
			# ignore genes with only 1 transcript
			my $cnt_transcripts = scalar keys %{$counts->{GENES}{$gene_id}};
			next unless $cnt_transcripts > 1;
			++$counts->{COUNT_GENES_MULTIPLE_TRANSCRIPTS};
			$counts->{COUNT_MULTIPLE_TRANSCRIPTS} += $cnt_transcripts;
	
			# save genes with multiple transcripts
			my $multi_trans_gene = \%{$counts->{GENES_WITH_MULTIPLE_TRANSCRIPTS}{$gene_id}};
			my @ensp = keys %{$counts->{GENES}{$gene_id}};
			@{$multi_trans_gene}{@ensp} = ();
		}


		if (exists $counts->{VIA_PARALOGS})
		{
			$counts->{COUNT_TRANSCRIPTS_VIA_PARALOGS} = scalar keys %{$counts->{VIA_PARALOGS}};
			{
				my %genes = map {$map_ensp_to_ensg->{$_} => $_} 
												keys %{$counts->{VIA_PARALOGS}};
				$counts->{COUNT_GENES_VIA_PARALOGS} = scalar keys %genes;
			}								
			$counts->{COUNT_TRANSCRIPTS_DISCARDED_OVERLAPPING} = scalar keys %{$counts->{DISCARDED_OVERLAPPING}};
			{
				my %genes = map {$map_ensp_to_ensg->{$_} => $_} 
												keys %{$counts->{DISCARDED_OVERLAPPING}};
				$counts->{COUNT_GENES_DISCARDED_OVERLAPPING} = scalar keys %genes;
			}								
			$counts->{COUNT_TRANSCRIPTS_DISCARDED_ORPHANS} = scalar keys %{$counts->{DISCARDED_ORPHANS}};
			{
				my %genes = map {$map_ensp_to_ensg->{$_} => $_} 
												keys %{$counts->{DISCARDED_ORPHANS}};
				$counts->{COUNT_GENES_DISCARDED_ORPHANS} = scalar keys %genes;
			}								
		}
	}
}


sub print_transcript_gene_counts(*\%\%)
{
	my ($FH, $counts_all_species, $map_ensp_to_ensg) = @_;
	print $FH "\t\t", $counts_all_species->{COUNT_SETS}, "\tortholog clades from \n";
	print $FH "\t\t", $counts_all_species->{COUNT_ORIG_CLADE}, "\toriginal transcript clades\n";
	my @species_names = ($ortho_name1, $ortho_name2);
	for my $species(0,1)
	{
		my $counts = $counts_all_species->{$species};
		print $FH "\t$species_names[$species]\n";


		#
		# genes and paralog before mergeing / disambiguating 
		#
		print $FH "\t\t", $counts->{COUNT_ORIG_TRANSCRIPTS}, 
					"\toriginal transcripts corresponding to ";
		print $FH $counts->{COUNT_ORIG_GENES}, 
					" genes\n";


		#
		# genes and paralog after mergeing / disambiguating 
		#
		if ($counts->{COUNT_TRANSCRIPTS})
		{
			print $FH "\t\t", $counts->{COUNT_TRANSCRIPTS}, 
						"\ttranscripts corresponding to ";
			print $FH $counts->{COUNT_GENES}, 
						" original genes\n";
		}


		#
		# resolve within clade specialties
		#
		if (exists $counts->{COUNT_TRANSCRIPTS_VIA_PARALOGS})
		{
			print $FH "\t\t", $counts->{COUNT_TRANSCRIPTS_VIA_PARALOGS}, 
						"\ttranscripts whose paralog has an ortholog, from ";
			print $FH $counts->{COUNT_GENES_VIA_PARALOGS}, 
						" original genes\n";
			print $FH "\t\t", $counts->{COUNT_TRANSCRIPTS_DISCARDED_OVERLAPPING}, 
						"\tdiscarded overlapping transcripts, from ";
			print $FH $counts->{COUNT_GENES_DISCARDED_OVERLAPPING}, 
						" original genes\n";
			print $FH "\t\t", $counts->{COUNT_TRANSCRIPTS_DISCARDED_ORPHANS}, 
						"\tdiscarded orphaned transcripts, from ";
			print $FH $counts->{COUNT_GENES_DISCARDED_ORPHANS}, 
						" original genes\n";
		}

		#
		# merged genes
		#
		if ($counts->{COUNT_GENES_MULTIPLE_TRANSCRIPTS})
		{
			print $FH "\t\t", $counts->{COUNT_GENES_MULTIPLE_TRANSCRIPTS},
						"\tmerged genes from $species_names[$species] split to give ";
			print $FH $counts->{COUNT_MULTIPLE_TRANSCRIPTS},
						" genes.\n";
			print $FH "\t\t", $counts->{COUNT_MULTIPLE_TRANSCRIPTS} - 
							$counts->{COUNT_GENES_MULTIPLE_TRANSCRIPTS},
						"\textra demerged genes.\n";
			print $FH "\t\t", $counts->{COUNT_EXTRA_TRANSCRIPTS},
						"\textra genes from $species_names[$species] altogether.\n";
		}
	
		#
		# orphaned genes
		#
		if ($counts->{COUNT_GENES})
		{
			print $FH "\t\t", $counts->{COUNT_ORIG_GENES} - 
								$counts->{COUNT_GENES},
						"\torphaned genes ($counts->{COUNT_ORIG_GENES} - $counts->{COUNT_GENES}).\n";
		}
	}

}

#_________________________________________________________________________________________
#
#	print_merged_paralogs
#
#		print out paralogs representing merged genes which have been disentangled
#
#_________________________________________________________________________________________
# FUNC_18
sub print_merged_paralogs(\%$)
{
	my ($counts_all_species, $file_name) = @_;

	# print species paralogs together
	open_or_die *MERGED, ">$file_name";
	print STDERR "$file_name\n";
	for my $species(0,1)
	{
		my $genes_with_multi = $counts_all_species->{$species}{GENES_WITH_MULTIPLE_TRANSCRIPTS};
		for my $gene_id(keys %$genes_with_multi)
		{
			my @transcripts = keys %{$genes_with_multi->{$gene_id}};
			print MERGED $gene_id, "\t";
			print MERGED '[', join(",", @transcripts), "]\n";
		}
	}
}


sub debug_dump_clades(\@$)
{
	my ($clades, $state) = @_;

	# no buffering of STDOUT
	select((select(STDOUT), $| = 1)[0]);
	for (0..$#$clades)
	{
		my $all_ensp = dump($clades->[$_]);
		$all_ensp =~ tr/\n//d;
		print $state." clade_id = $_, $all_ensp\n";
	}
}





#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Ensp to ensg from panda
#
#
print STDERR "\t",  "Get ensp and ensg from panda\n";
my %map_ensp_to_ensg;
# FUNC_1
get_ensp_to_ensg($dbh, %map_ensp_to_ensg);




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Read in prot coding loci from panda
#
#
print STDERR "\t",  "Get ensp protein coding loci\n";
my %map_ensp_to_loci;
# FUNC_2
get_prot_coding_loci($dbh, $ortho_name1, $ortho_name2, %map_ensp_to_loci);



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#
#	initialize paralog results
#


use constant INITIAL		=> 0;
use constant RESOLVED		=> 1;
use constant SINGLE_CHOICE	=> 2;
use constant UNAMBIGUOUS	=> 3;
use constant AMBIGUOUS		=> 4;
use constant TOTAL			=> 5;
my @names = qw(initial resolved single_choice unambiguous ambiguous total);
my @ortholog_results;
push(@ortholog_results, {}) for (0..5);
for my $i(0..$#ortholog_results)
{
	# FUNC_3
	init_paralog_result_struct(%{$ortholog_results[$i]}, $names[$i]);
}





#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Read in dS pairs
#
#
print STDERR "\tRead in dS pairs\n";
my %seq_pairs_dS;
read_ds_pairs($ds_pairwise_file_name, %seq_pairs_dS, undef);







#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Read in ortholog clades
#
#
my @all_clades;
print STDERR "\tRead in ortholog clades from $ortholog_clades_file_name...\n";
# FUNC_4
read_in_ortholog_clades($ortholog_clades_file_name, @all_clades);





#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Split clades by species
#
#
# FUNC_5
print STDERR "\t",  "Split transcripts by species\n";
split_transcripts_by_species(@all_clades);




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Save what we have so far
#
# FUNC_6
print STDERR "\t",  "Record initial transcripts and gene set\n";
record_unclustered_clades(@all_clades, %{$ortholog_results[INITIAL]});

###DEBUGGG
#print "Unclustered clades\n", dump @ortholog_results, "\n\n";


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Resolve ensp per orthologue clade of transcripts
#		each orthologue clade must not contain overlapping transcripts from multiple
#			genes
#		also calculate mean ortholog ds and splits transcripts by species
#
#
my %resolve_clades_results;
init_paralog_result_struct(%resolve_clades_results, "Resolve per clade results");
#	DISCARDED_OVERLAPPING
#	DISCARDED_ORPHANS
#	VIA_PARALOGS
# FUNC_7
print STDERR "\t",  "Resolve overlapping transcripts within single clades\n";
resolve_overlapping_transcripts_within_all_clades(	@all_clades,
													%map_ensp_to_ensg, 
													%map_ensp_to_loci, 
													%seq_pairs_dS,
													%{$ortholog_results[RESOLVED]});


###DEBUGGG
#debug_dump_clades(@all_clades, "RESOLVED");

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	calculate mean dS for each clade
#
#
# FUNC_8
my @mean_ortholog_ds;
print STDERR "\t",  "Calculate mean dS to orthologs for all clades\n";
calculate_mean_ortholog_dS_for_all_clades(@all_clades, %seq_pairs_dS, @mean_ortholog_ds);




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	in case some ensp still belong to multiple clades
#		pick clade with lowest dS
#
#
# FUNC_9
print STDERR "\t",  "Assign each transcript to the clade with the lowest dS\n";
assign_ensp_to_clade_with_lowest_ds(@mean_ortholog_ds, @all_clades);

# some clades have changed: recalculate mean ortholog dS
@mean_ortholog_ds = ();
print STDERR "\t",  "Re-calculate mean dS to orthologs for all clades\n";
calculate_mean_ortholog_dS_for_all_clades(@all_clades, %seq_pairs_dS, @mean_ortholog_ds);



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	discard clades without orthologous matches
#
{
	my @only_orthologous_clades;
	my @only_orthologous_clades_ds;
	for my $clade_id(0.. $#all_clades)
	{
		next if $mean_ortholog_ds[$clade_id] >= 10.0;
		push(@only_orthologous_clades,  $all_clades[$clade_id]);
		push(@only_orthologous_clades_ds, $mean_ortholog_ds[$clade_id]);
	}
	@all_clades = @only_orthologous_clades;
	@mean_ortholog_ds = @only_orthologous_clades_ds;
}



###DEBUGGG
#debug_dump_clades(@all_clades, "LOWEST");
#for my $clade_id(0.. $#all_clades)
#{
#	print "LOWEST clade_id = $clade_id, ds = $mean_ortholog_ds[$clade_id]\n";
#}




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Save what we have so far
#
print STDERR "\t",  "Record transcripts and gene set after removal of overlapping\n";
record_unclustered_clades(@all_clades, %{$ortholog_results[RESOLVED]});


	
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	translate ensp to ensg
#
#
my @all_clade_genes;
my %ensg_to_clade_id;
for my $clade_id (0 .. $#all_clades)
{
	my ($ensp1, $ensp2) = @{$all_clades[$clade_id]};
	my %ensg1 = map {$map_ensp_to_ensg{$_} => undef} @$ensp1;
	my %ensg2 = map {$map_ensp_to_ensg{$_} => undef} @$ensp2;
	push(@all_clade_genes, [ [keys %ensg1], [keys %ensg2] ]);
	for (keys %ensg1, keys %ensg2)
	{
		push(@{$ensg_to_clade_id{$_}},  $clade_id);
	}
}



##DEBUGGG
#select((select(STDOUT), $| = 1)[0]);
#for (0..$#all_clade_genes)
#{
#	my $all_ensp = dump($all_clade_genes[$_]);
#	$all_ensp =~ tr/\n//d;
#	print "GENES clade_id = $_, $all_ensp\n";
#}


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	cluster ortholog clades which have shared genes together
#
#
print STDERR "\tCluster ortholog clades which have shared genes together\n";
my (@cross_linked_clades, @single_choice_clades);
# FUNC_10
print STDERR "\t",  "Cluster clades by shared genes\n";
cluster_ortholog_clades(@cross_linked_clades, @single_choice_clades, 
						scalar @all_clades, %ensg_to_clade_id);



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	Divide cross_linked_clades into unambiguous and ambiguous clades
#
#
my @unambiguous_clades;
my @ambiguous_clades;
for my $clade_set (@cross_linked_clades)
{
	# FUNC_11
	if (ortholog_clades_in_cluster_are_identical(@$clade_set, @all_clade_genes))
	{
		push(@unambiguous_clades, $clade_set);
	}
	else
	{
		push(@ambiguous_clades, $clade_set);
	}
}










#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	print out orthologs with no multiple transcript options
#
#
# FUNC_12
{
	print STDERR "\t",  "Print out clades without transcripts only present in single clades\n";
	resolve_singly_spliced_transcripts(	@single_choice_clades, 
										@mean_ortholog_ds,
										@all_clades,
										%map_ensp_to_ensg,
										%seq_pairs_dS,
										%{$ortholog_results[SINGLE_CHOICE]});
	$ortholog_results[SINGLE_CHOICE]{COUNT_ORIG_CLADE} = 
	$ortholog_results[SINGLE_CHOICE]{COUNT_SETS}		= scalar @single_choice_clades;
	print STDERR "\tSingle transcript clades printed\n";
}








#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	print out unambiguous clades
#
#
{
	#
	#	handle merge genes within the same clade for "unambiguous" clades
	#
	print STDERR "\tPrint out unambiguous clades.\n";
	open_or_die *FH_EASY_RESOLVED, ">$unambiguous_file_name";
	# FUNC_13
	count_clades_in_sets(%{$ortholog_results[UNAMBIGUOUS]}, @unambiguous_clades);
	for my $clade_set(@unambiguous_clades)
	{
		# FUNC_14
		resolve_unambiguous_clades(	$clade_set,  
									@mean_ortholog_ds, 
									@all_clades, 
									%map_ensp_to_ensg,
									*FH_EASY_RESOLVED,
									%seq_pairs_dS,
									%{$ortholog_results[UNAMBIGUOUS]});

	}
	print STDERR "\tUnambiguous transcript clades printed\n";
}







#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#	print out ambiguous clades
#
#
{
	#
	#	handle merge genes within the same clade for "ambiguous" clades
	#
	print STDERR "\tPrint out ambiguous clades.\n";
	open_or_die *FH_HARD_RESOLVED, ">$ambiguous_file_name";
	
	count_clades_in_sets(%{$ortholog_results[AMBIGUOUS]}, @ambiguous_clades);
	for my $clade_set(@ambiguous_clades)
	{
		# FUNC_15
		resolve_ambiguous_clades(	$clade_set,  
									@mean_ortholog_ds, 
									@all_clades,
									*FH_HARD_RESOLVED,
									%map_ensp_to_ensg,
									%map_ensp_to_loci,
									%seq_pairs_dS,
									%{$ortholog_results[AMBIGUOUS]});

	}
	print STDERR "\tAmbiguous transcript clades printed\n";
}



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#
#	print SUMMARY
#

#	aggregated results
for my $ortholog_results(@ortholog_results[SINGLE_CHOICE, UNAMBIGUOUS, AMBIGUOUS])
{
	# FUNC_16
	merge_paralog_result_struct(%{$ortholog_results[TOTAL]}, %$ortholog_results);
}

###DEBUGGG
#print dump @ortholog_results, "\n\n";


print STDERR "\tPrint out summary.\n";
{
	open_or_die (*SUMMARY, ">$summary_file_name");
	for my $ortholog_results(@ortholog_results)
	{
		my $name = $ortholog_results->{NAME};
		print SUMMARY "\n$name Orthologs Clades\n";
		# FUNC_17
		finish_transcript_gene_counts(%$ortholog_results, %map_ensp_to_ensg);
		print_transcript_gene_counts(*SUMMARY, %$ortholog_results, %map_ensp_to_ensg);
		# FUNC_18
		print_merged_paralogs(%$ortholog_results, "$dir_orthologs_output/resolve_orthologs/$name.merged.paralogs");
	}
}
system ("cat $summary_file_name");






#
#   print list of proteins with no loci
#
if (scalar keys %ensp_sans_loci)
{
    print STDERR "\t", scalar keys %ensp_sans_loci, " sequences without loci data\n";
    print STDERR "\tSee $sans_loci_error_file_name\n";
    open_or_die *ERR_OVERLAP, ">$sans_loci_error_file_name";
    print ERR_OVERLAP join ("\n", keys %ensp_sans_loci);
}

log_pipeline_stage($dbh, $curr_protocol_id, 31, 'Resolve phyop orthologs', $start_time);
$dbh->disconnect();



#TODO clusters of empty ignored
#TODO no ortholog dS ignored

